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ПРИ ЧИСЛЕННОМ ИНТЕГРИРОВАНИИ ЗАДАЧИ 
ИЗОЭЛЕКТРИЧЕСКОГО ФОКУСИРОВАНИЯ 


Л.В. САХАРОВА 
(Филиал Морской государственной академии имени адмирала Ф.Ф. Ушакова в г. Ростове-на-Дону) 


Проведен математический анализ проблемы возникновения так называемых «негауссовских» режимов при 
численном решении интег/ ро-дифференциальной задачи изоэлектрического фокусирования. В результате 
аналитического преобразования задачи, составления оптимизационных алгоритмов и численного 
асимптотического тестирования установлено, что «негауссовские» режимы являются СВОЙСТВОМ ИСХОДНОЙ 
математической задачи, а не результатом накопления вычислительной пог, решности. 

Ключевые слова: математическое моделирование, краевая задача, «негауссовские» режимы, гауссовское 
распределение, оптимизационный алг Ъоритм. 


Введение. Гауссовское (нормальное) распределение является одним из важнейших и наиболее 
часто встречающихся распределений в естественных науках. Теория гауссовского распределения 
лежит в основе математической теории изоэлектрического фокусирования (ИЭФ) — эффективного 
и универсального метода фракционирования и анализа белков [1]. Практика подтверждает, что 
гауссовское распределение применимо к широкому классу так называемых амфолитов-носителей, 
т. е. амфотерных кислот с высокой буферной емкостью. Распределение концентраций компонент 
смеси имеет гауссовский вид: 


С= С, ехр(-рЕх? /25), 
где С - концентрация; Е - напряженность поля; Р - коэффициент диффузии; р - градиент 


ы Чи 
электрофоретической подвижности амфолита, р=-—. 


фх 
Гауссовское распределение концентраций амфолитов многими зарубежными авторами 
было получено при компьютерном моделировании ИЭФ [2 - 4], а также искажение гауссовского 
распределения (рис. 1) [5 - 7]. 
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Рис. 1. Негауссовские (а) и классические гауссовские режимы (6) [6] 
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При высоких плотностях электрического тока распределение концентраций на графиках 
приобретает «платообразную форму», которая существенно отличается от классического гауссов- 
ского (а также распределений Хи-квадрат, Стьюдента, Коши, Вигнера и пр.). При высоких плотно- 
стях тока соответствующая начально-краевая задача приобретает ряд особенностей, затрудняю- 
щих решение и приводящих к неконтролируемому накоплению вычислительной погрешности. По- 
этому возникает закономерный математический вопрос (оставшийся за рамками работ [5 - 7], 
являющихся прикладными электрохимическими исследованиями): не обусловлены ли наблюдае- 
мые «гипергауссовские» режимы результатом систематического накопления погрешности? Кор- 
ректен ли математически вывод о существовании «негауссовских» режимов ИЭФ? Ответить на 
этот вопрос помогло представленное в работе исследование математической модели ИЭФ анали- 
тическим и численным методами. 

1. Физическая и математическая постановка задачи. Дано: водный раствор № амфолитов 
помещается в электрофоретическую камеру (ЭК), представляющую собой цилиндр длиной / и 


радиусом »х. Для каждого амфолита известны его коэффициенты миграции и, , константы диссо- 
циации реакций К“ и К®, атакже общие количества т, , К=1,...,М. Температура Т’внутри ЭК 


считается постоянной. Под действием постоянного тока плотности „ЛГ в ЭК формируется распре- 
деление концентраций амфолитов, ведущее к стационарному распределению концентрации ионов 
водорода рН’. В данной модели рассматривается продольное осевое сечение ЭК, представляю- 


щее собой прямоугольный цилиндр длиной / и шириной 2х (рис. 2). 




















Рис. 2. Распределение концентрации трех амфолитов при ИЭФ: а) профили концентраций и графики рН 
и внутренней проводимости ЭК; 6) продольное (осевое) сечение ЭК со стационарным распределением 


Предполагается, что в водном растворе диссоциация „ -го амфолита протекает по схемам: 
МН+ВСООН <> МН,ВСООН +Н*, МН,ВСООН <> МН,ВСОО' +Н*, 
где  МНВСООН, МН,ВСОО`, МН.,ВСООН - положительный, отрицательный и нейтральный 
ионы амфолита, молярные концентрации которого обозначены &, &, &. 
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Общая или так называемая аналитическая концентрация амфолита определяется форму- 
лой &, = +8 +“. Помимо приведенных реакций диссоциации амфолитов в водном растворе 
следует учитывать реакцию автопротолиза воды: Н.О > ОН`+Н“*, отсюда следует, что: 


2 
ОН 5. где Н - концентрация ионов водорода, ОН `- концентрация гидроксил-ионов и К. = 


константа автодиссоциации воды, Е =10-"“. В сделанных предположениях (при отсутствии пото- 


НА 
ка вещества через границы ЭК) основное уравнение теории переноса И. = 0 приобретает 
р 


упрощенную форму: 1 =0. Поэтому стационарное распределение концентраций &, (х) и распре- 
деление рН(х) описывается следующей одномерной интегро-дифференциальной задачей (иско- 
мые функции не зависят от ординаты, и картина неизменна в любом осевом сечении ЭК): 





_ Ч +6, (9-0 )Е=0, К=1,2,..,М; (1) 
ах 
у а ; К К АН 
7-52 к (© -а^ 5, + (о на! )ьЕ |2, —— инНЕ + 
а(ОН 
+Бон О) } НН рОНЕ; (2) 
ах 

М 
Уи - Е, +Н-ОН =0; (3) 

К=1 

1 
ле? [6 сдак=т,,  К=12,..М, (4) 


Здесь = - стандартный электрохимический параметр, == АТ/ЁР, где величины К, Т и Е- уни- 
версальная газовая постоянная, температура и число Фарадея соответственно; и„, но» -— под- 


вижности ионов водорода и гидроксил ионов; О,, Р‚„, О», - коэффициенты диффузии ионов, 
р; =4, , &- степень диссоциации амфолита. 


Уравнение массопереноса получено на основании уравнения потока амфолита (1). Соот- 
ношение (2) представляет собой обобщенный с учетом диффузии закон Ома (плотность тока яв- 
ляется суммой плотностей токов всех ионов, включая ион водорода и ион гидроксила), (3) есть 
уравнение электронейтральности. Наконец, интегральное условие (4) соответствует закону со- 
хранения массы вещества, в соответствии с которым суммарное количество всех трех форм ам- 
фолита неизменно и равно т,. 


Основная математическая трудность численного интегрирования системы (1) - (4) заклю- 
чается в том, что: 

а) при решении системы дифференциальных уравнений (1) относительно концентраций 
необходимо определять величину Н из алгебраического уравнения (3); 

6) вместо обычных краевых условий приходится использовать интегральные условия (4). 

Тем не менее, для низких и средних значений „Л задача решается стандартными асимпто- 
тическими и численными методами и дает классические гауссовские распределения концентраций 
[8] - [10]. 

Для больших значений „Л, как следует из (1), перед функциями &, появляется большой 


параметр Л( (== 0,0257 ). Это приводит к дополнительным трудностям: 
в) слабые изменения &, вызывают существенные изменения производных 4&, /4х , что 


может привести к неконтролируемому накоплению вычислительной погрешности; 

Г) в «гипергауссовских» режимах (рис. 3, д‚е) в области «плато» (а также за пределами 
области фокусирования амфолита) производные расчетных функций близки к нулю и вероятно 
«зацикливание» метода Рунге-Кутта, сопровождаемое выходом на неверное решение-константу; 


7 


Физико-математические науки 








д) в остальных точках производные расчетных функций, наоборот, стремятся к бесконеч- 
ности, что может вызвать резкий скачок решения вплоть до выхода его на отрицательные (ли- 
шенные физического смысла) решения; 

е) за пределами области фокусирования амфолита его концентрация стремится к нулю, 
что может привести к обнулению искомых функций при недостаточной точности вычислений. 

Таким образом, предварительный анализ задачи и имеющихся результатов численного 
решения задачи математического моделирования ИЭФ в «негауссовских» режимах позволяет сде- 
лать вывод: для задачи в стандартной формулировке достаточно велик риск неконтролируемого 
накопления погрешности. Как следствие возникает вопрос: существуют ли, с математической точ- 
ки зрения, «негауссовские» решения задачи (1) - (4) или они являются результатом систематиче- 
ского накопления погрешности? Чтобы ответить на этот вопрос, необходимо объемное аналити- 
ческое исследование. 

2. Преобразование системы. Теорема 1. Интегро-дифференциальная задача (1) - (4) относи- 
тельно №М-+2 неизвестных функций Н, Е, &,(х), К=1,2,...,№М может быть сведена к краевой 


задаче относительно 2№ неизвестных функций а,(х), п, (х), К=1,2,...,М: 
= Ча 1 = Ф, (у) У 





фа, (у) о’ в 
т (у), (6) 

ах 
_х не ИИ о еы й 
б Фр (ФС) А у+исй(у-ч,), (7) 
Ф, (у) = 5, + сй(у — у, Г (8) 
у= и | + У очки.) = 5 + У, есчо) я (9) 
п,(0)=0, (10) 
т =, &=1,2,...М. (11) 

дг 

При этом 

&,(х)=2Ка,(х)ф, (у), (12) 
Н=К,ехр(у). (13) 


Доказательство. На первом этапе доказательства в рассмотрение вводится новая функция чу 
(определена уравнением (13)), а для упрощения уравнений - новые параметры: 





1 (кКФк® 

=5М |, 14 

ЗК” 
5, в- Ко р (15) 
че, (16) 

2 Ин 
И-УИнИон ' (17) 

В новых обозначениях функции, входящие в (1) - (3), приобретают форму: 
51 (ум, сп(у-,) 








р +, = 6, = 
6, +ей(у-\,) б, + ей(у-у,) 
Кроме того, используются новые функции и новая плотность тока: 

век", ЛЕК, Л, 


р р 
о -а=е, = 
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В результате, опуская символ ипем’, систему (1) - (3) можно переписать в форме, 
обладающей существенным преимуществом: более компактный вид плюс отсутствие малого 
параметра А, дающего разброс слагаемых по степеням малости: 





—5 Ч, +, е Е =0, (18) 
ах 
ый а 
= А [5 6, ) +0, Е д (-=Уу Е Е)нсй(у =) (19) 
= 
М 
Ув, Е: +5йу=0. (20) 


&=1 
На втором этапе доказательства для упрощения системы вводятся новые функции 
(определены как (8)) : 








а а? 
та (=, ом = =), 
да Фи (у) _ ФЕ Си) (21) 


р [1 ; 
Ф, (у) Фу) 

Важно, что функции &,(х) могут быть найдены в следующей форме: &, (х) = а. ф,(у), где 
а, (х) — новая неизвестная функция. В новых переменных система (18) - (21) будет иметь вид: 


Ча, 1 _Ф ` „_ ($ 
Е, о = а (фи) + исй(у фм), 
О ор .(Ф _ )+ исй(у — уо) 


(определены как (5) и (7)) 


Ха Фу +эйу=0. (22). 
Е=1 
Система из (5), (7), (22) значительно проще, чем (18), (19), (20) , так как из нее исключена 
неизвестная функция Е. По сути, уравнения (5) и (7) могут быть записаны одним 
дифференциальным уравнением. 
На третьем этапе (22) получает дополнительное преимущество: оно разрешается 
относительно функции зу. Предположим, что 


Ру) = У зи ур) +5йч , (23) 


тогда Р(Су)=Р”Су). Решением последнего уравнения, очевидно, является функция 
Е (у) = Азйу + Взрлу . Из (23) следует 


Е(0)= А=-У а, 5йу,, Е'(0)=В=-У а, сйу, +1. 
= = 


(22) эквивалентно уравнению Р(у)=0 . Следовательно, лу =-— А/В , а значит, 


1 Ш 1 Ш 
у= 5. | +», очки.) Е 5 +У а, есчо) | 
Е=1 1 


К= 
(обозначено как (9)). 
На четвертом этапе интегральные условия (4) трансформируются в обычные краевые 
условия посредством введения новых функций, при этом ликвидируется главная вычислительная 
трудность задачи 


п‚(®) = [`а,ф, (уда, (24) 
удовлетворяющие следующей краевой задаче ( обозначена как (6), (10), (11)): 
ап, (х т 
419 (у), п (0)=0, и) =. 
фх ди 


Теорема доказана. В результате избавляемся от трудностей численного интегрирования задачи 
(а) и (6), указанных в п. 1. 
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3. Алгоритмы численного решения задачи. Чтобы избежать указанных в п. 1, в -е проблем, 
производим дополнительную замену переменной, составление оптимизационных алгоритмов с 
контролем выхода на неверные решения. 


Неизвестная функция с, может быть представлена в виде экспоненты 
с, =Ь, ехр(И8Е, (х)), Е=1,2....М, (25) 
где Б, — постоянный параметр. Экспоненциальная форма решения «отсекает» выход на лишен- 
ные физического смысла отрицательные решения. Сомножитель 1/= обеспечивает «сглажива- 
ние» расчетной функции: малым изменениям функции А, (х) соответствуют существенные изме- 
нения функции с,. Для неизвестных функций 2, (х), &=1,2,...,№М может быть получена краевая 
задача с условиями (10), (11): 


АЕ, _ Ф(у) Л (6) 
а  ф,(у)о 
ие = ф, (У)Б, ехр [1 69) ; (27) 
-Х 5 
с= > (оу) о, ехр [1 6 Е 28) 


1 п 1 1 Ш 1 
\у — —т 1 + У ехр [к (х) т м, — —т 1 тр ХА оч (х) га м, Й (29) 
2 = 5 2 Га 5 
где параметры и,, и’, \,, у, определяются формулами (14) - (17), а функции &, (х) и Н - из 


уравнений (8) ‚ (12), (13), (25). 

Для решения задачи (26) — (29), (10), (11) разрабатываются два оптимизационных 
алгоритма, обеспечивающих необходимую точность решения. 

Первый оптимизационный алгоритм основан на модифицированных методах Рунге-Кутта и 
Ньютона и обеспечивает эффективное решение начально-краевой задачи при заданной плотно- 
сти тока „Л. Он состоит из трех этапов. На первом этапе задаются начальные приближения: 


Е, (0)=х, , К=1, 2... М, (30) 
а также решаются задачи (26) - (29), (10), (30) на отрезке [0,1] методом Рунге-Кутта (но полу- 
ченные при этом значения и, (Г) не удовлетворяют условиям (11)). На втором этапе модифициро- 
ванным методом Ньютона осуществляется поиск точных значений х, . На третьем этапе выполня- 
ется решение задачи (26) - (29), (10), (30) для точных значений х, ‚ полученных на предыдущем 
этапе; одновременно вычисляется &,, рН =-1№Н и о, а также строятся графики соответст- 
вующих функций. 

Сложность работы по первому оптимизационному алгоритму состоит в выборе начальных 
значений 2; (0) = х,, достаточно близких к точным. Численный эксперимент показал, что произ- 
вольный выбор этих значений (для больших значений ./) часто приводит к расхождению итера- 
ций метода Ньютона либо получению неверных решений. В то же время для малых величин .Л 
начальные приближения 2, (0)=х, легко определяются из асимптотических формул и всегда 
приводят к верному результату. Поэтому для исследования задачи во всем диапазоне „/ исполь- 
зуется так называемый метод движения по параметру, когда величина „Г возрастает с перемен- 
ными контролируемыми шагами, а за начальные приближения принимаются точные значения, 
полученные на предыдущем шаге 2, (0) =х,. 


Второй оптимизационный алгоритм, реализующий метод движения по параметру, позво- 
ляет выполнять расчеты в широком диапазоне плотностей тока, выбирать оптимальный шаг дви- 
жения по плотности тока, а также контролировать процесс накопления вычислительной погреш- 
ности, отсекая при этом возможные выходы на неверные решения (возникающие в случае выбора 
слишком большого шага). Он также состоит из трех этапов. На первом этапе задаются параметры 
ячейки и характеристики амфолитов, на основании заданных величин рассчитываются постоян- 
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ные значения ф,, 6,, Ф,, и и определяются начальные приближения х, на основе асимптоти- 
ческих формул. 
На втором этапе к начальной плотности тока „Л, применяется Алгоритм 1, строятся гра- 


й 


фики &,, рНи с. Наконец, на третьем этапе выполняется циклическое движение по плотности 
тока „Л =/Л-+ АЛ с применением на каждом шаге Алгоритма 1. При вычислениях осуществляется 
постоянный контроль условия: тах &, (х)>&, А=1,2,...,№, невыполнение которого означает 


выход на неверное решение, а также необходимость уменьшить шаг по плотности тока и, 
вернувшись на шаг назад, выполнить повторный расчет. 

Алгоритмы реализуются на языке программирования ТигБо Разса! с использованием 
стандартного модуля Сгарп. 
4. Тестирование модели асимптотическими методами. Модель обнаруживает качественное 
соответствие с общепринятыми математическими моделями ИЭФ при низких и средних плотно- 
стях тока. Получаемые посредством нее профили концентраций амфолитов имеют вид стандарт- 
ных гауссовских распределений (рис. 3). При высоких плотностях тока, как и в [6; 7], системати- 
чески наблюдаются «негауссовские» «платообразные» режимы. Вначале максимумы на профилях 
трансформируются в «плато», а затем сами профили приобретают вид прямоугольников, вплот- 
ную примыкающих друг к другу. 


и и 
ммоль/л ммоль/л 








Хх, ДМ 





а) 


и 
ммоль/л ммоль/л 








Рис. 3. ИЭФ системы из пяти стандартных амфолитов с рН>7: 
Н&-—Н5, Ны-ОСр, Н5, В- Аа -Нь , Туг- Аго 
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Однако все принятые меры к минимизации вычислительной погрешности еще не дают ос- 
нования с полной уверенностью говорить, что в данном случае удалось избежать систематиче- 
ских ошибок и наблюдаемые «негауссовские» режимы существуют в действительности. 

Поэтому для проверки адекватности построенной расчетной модели был использован 
асимптотический метод тестирования. Он базируется на очевидном утверждении. 

Утверждение. Если при рассматриваемой плотности тока „Л система (26) - (29) имеет 


«платообразное» («негауссовское») решение, то на отрезке м соответствующем области 
фокусирования А-го амфолита, его концентрация &, и функция у постоянны (рис. 3): 
г хЕ[х х ] 
— А > х-12 7” и 
6, = ‚ \-ЕУ,, хе [Еж] . 
0,хя [5% ] 


Доказана теорема 2. Если при рассматриваемой плотности тока „Л система (26) - (29) 
имеет «платообразное» («негауссовское») решение, то значения функций Е, (х) при х=0 опре- 
деляются асимптотической формулой: 








К 
Е, (0) = ша -У/У,Ф,(у,а?), (33) 
11 
где 
(у) 1 
Ф-Т, (34) 
ф, (у) о 
п ее (35) 
й, (9. +1) 2К г 


вт [У] (36) 


5. Результаты расчетов и их интерпретация. Расчеты проводились в предположениях: длина 
ЭК /=2 дм; радиус ЭК „=0,2 дм; Т=298 К. Плотность тока измеряется в А/дм?. Рассмотрим 


систему из пяти стандартных амфолитов с рН>7:Н5-Ны, Н&-С, НБ, В-Аа-Ны, 
Туг- Аго (табл. 1) [1]. Исходные количества амфолитов одинаковы, М, =0,1 моль. Имеет место 
существенная неравномерность распределения амфолитов по константам диссоциации К“, К“® 




















К®+рк® 
и изоэлектрическим точкам: р/ В (табл. 1). 
Таблица 1 
Характеристики амфолитов 
№ п/п Амфолит рК® рК® р ДрКк Подвижность, у-10^* 
1 Нё-Нб 6,80 7,80 7,30 1,00 1,49 
2 Нё-Су 6,27 8,57 7,42 2,30 2,40 
Е Нр 6,00 9,17 7,59 3,17 2,85 
4 В- А/а- Н5 6,83 9,51 8,17 2,68 2,30 
5 Туг - Аге 7,55 9,80 8,68 2,25 1,58 





























При плотности тока .Л =0,0005 расслоение незначительно, и ни один из амфолитов не 
достиг еще своего изоэлектрического состояния (рис. 3, а). При плотности тока ./ =0,002 на про- 


филе В - Аа — Ни отчетливо проявляется максимум, указывающий на достижение амфолитом 
изоэлектрической точки (рис. 3, 6), прослеживается появление максимума на профиле Н&. При 
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/=0,007 (рис. 3, в) изоэлектрического состояния достигают все пять амфолитов, причем про- 
филь В - Аа - Нк имеет вид стандартного гауссовского распределения, в то время как профили 
Нб — б/у и Нб асимметричны; график рН монотонно возрастает. При Л =0,027 (рис. 3, Г) на про- 
филе Туг- Ага появляется «плато», он утрачивает сходство с гауссовским распределением, а 
значит, система входит в «негауссовский» режим, характеризующийся «негауссовскими» распре- 
делениями концентраций. При .7=0,235 (рис. 3, д) такие «плато» видны уже на всех профилях; 
амфолиты расслаиваются на полосы одинаковой ширины; профили имеют вид правильных пря- 
моугольников для Н, В — Аа — Н&, Туг — Ато. Наконец, при Л =1,323 (рис. 3, е) профили всех 
амфолитов имеют вид прямоугольников, внутри которых рН и с постоянны. 

Для данной системы проводилось асимптотическое тестирование по формуле (20) 
(табл. 2), сравнивались расчетные значения 2, (х), А =1,2,..., № (т. е. полученные в результате 
работы основной программы) и асимптотические, рассчитанные в предположении, что имеет ме- 
сто «платообразный» («негауссовский») режим. 

Для Л =1,307 расчетные и асимптотические значения имеют расхождение в пятом знаке 
после запятой (табл. 2). Это указывает на высокую точность вычислений в «негауссовских» 
режимах. Поскольку вычисление расчетных величин А, (0) не зависит от вычисления 


асимптотических величин Е, (0) ‚ то можно сделать основной ВЫВОД: феномен «негауссовских» 


режимов является свойством математической (и соответствующей физической) задачи, а не 
результатом накопления вычислительной погрешности. 



































Таблица 2 
Расчетные и асимптотические значения А, (х) 
Величины Е(0) Е, (0) В(0) Е (0) Е. (0) 

= 0,075 

Расчетные —0,0066864 —0,0091958 —0,2125067 —2,9209963 —10,3920329 

Асимптотические —0,0072057 —0,1321356 —0,2679040 —3,3555129 —11,864537 
Л=0,347 

Расчетные —0,0066864 —0,2951293 —0,7671985 —13,3475905 —47,9548238 

Асимптотические —0,0067035 —0,3243769 —0,7910892 —13,187091 —7,272132 
Л =1,307 

Расчетные —0,0067712 —1,0472041 —2,7635703 —50,2134548 —180,505632 

Асимптотические —0,0067679 —1,0471994 —2,7635745 —50,213453 —180,505631 























Заключение. 1. Система дифференциальных уравнений с интегральными условиями (1) - (4) 
может быть преобразована к обычной краевой задаче, преимуществами которой являются: 

— замена интегральных условий краевыми посредством введения дополнительных 
неизвестных функций; 

— исключение из рассмотрения сложного алгебраического уравнения; 

— исключение из рассмотрения слагаемых с малым параметром Х.,. 


Указанные преимущества помогают свести к минимуму факторы, приводящие к 
накоплению погрешности при численном решении задачи. Оптимизирует поиск решения 
представление его в экспоненциальной форме, исключающей выход на неверные отрицательные 
решения, а также оптимизационные алгоритмы, основанные на методах Рунге-Кутта, 
модифицированном методе Ньютона, методе движения по параметру. 

2. Расчетная модель обнаруживает качественное соответствие с классическими 
гауссовскими моделями ИЭФ при низких и средних плотностях тока Л. Вычисления показывают, 
что при достижении некоторого критического значения „Л профили концентраций теряют 
сходство с гауссовским распределением. Вначале вершины кривых трансформируются в «плато»; 
затем профили превращаются в прямоугольники, примыкающие друг к другу. Графики рН и с 
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имеют ступенчатую форму, это «платообразные», «гипергауссовские» распределения, 
представляющие серьезный научный интерес. 

3. Тестирование расчетной модели на основе асимптотического метода показало высокую 
степень соответствия расчетных и асимптотических значений контрольной функции. 

Таким образом, можно констатировать существование у задачи математического 
моделирования ИЭФ «гипергауссовских» решений, в которые трансформируется классическое 
гауссовское распределение при достижении некоторой критической плотности тока. 
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